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Computational fluid dynamics (CFD) has the potential to improve the historical rocket 
injector design process by simulating the sensitivity of performance and injector-driven 
thermal environments to the details of the injector geometry and key operational 
parameters. Methodical verification and validation efforts on a range of coaxial injector 
elements have shown the current production CFD capability must be improved in order to 
quantitatively impact the injector design process. This paper documents the status of an 
effort to understand and compare the predictive capabilities and resource requirements of a 
range of CFD methodologies on a set of model problem injectors. Preliminary results from a 
steady Reynolds-Average Navier-Stokes (RANS), an unsteady Reynolds-Average Navier- 
Stokes (URANS) and three different Large Eddy Simulation (LES) techniques used to model 
a single element coaxial injector using gaseous oxygen and gaseous hydrogen propellants are 
presented. Initial observations are made comparing instantaneous results, corresponding 
time-averaged and steady-state solutions in the near-injector flow field. Significant 
differences in the flow fields exist, as expected, and are discussed. An important preliminary 
result is the identification of a fundamental mixing mechanism, accounted for by URANS 
and LES, but missing in the steady RANS methodology. Since propellant mixing is the core 
injector function, this mixing process may prove to have a profound effect on the ability to 
more correctly simulate injector performance and resulting thermal environments. Issues 
important to unifying the basis for future comparison such as solution initialization, 
required run time and grid resolution are addressed. 
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I. Introduction 

The rocket engine thrust chamber assembly (TCA) is the heart of a rocket engine in the sense that it must deliver 
a specified level of thrust, as required by the vehicle, to achieve its mission. The TCA includes the injector, 
combustion chamber, nozzle, and the associated manifolds that deliver the propellants and coolant. Successful and 
robust TCA designs must provide reliable ignition, achieve a high level of performance, produce manageable 
thermal environments, and exhibit highly stable combustion characteristics. The injector design has a significant 
impact on these competing requirements and thus represents one of the most critical components of the design cycle. 

Flow physics associated with rocket engine injectors are very complex. Phenomena associated with propellant 
injection are sensitive to the thermodynamic state of the fuel and oxidizer streams and include complex processes 
such as atomization, vaporization, turbulent mixing, and combustion. Propellants are often at cryogenic temperatures 
and very high pressures as they enter the injector. The resultant flow downstream of the injector is spatially 
evolving, highly unsteady, three-dimensional, and involves a wide range of scales in time and space. 

Historical injector design methodologies in the rocket propulsion community are characterized by relatively 
simple design tools 1 and large amounts of costly testing. The empirical design tools are either zero- or one- 
dimensional and do not adequately describe the injector geometry or complex flow phenomena that influence the 
operational performance of a TCA. The test, fail, fix cycle that is inherent in the historical injector development 
process is a direct consequence of these modeling inadequacies. As a consequence, there are many examples of 
combustion devices failures over the last 50 years. There were 18 major combustion devices failures in the Space 
Shuttle Main Engine development program alone. 2 

NASA’s Constellation Program has neither the schedule nor budget to accommodate the typical engine 
development program. This reality has led to a decision to modify existing engines, namely the J-2 and RS-68, to 
meet propulsion requirements for the ARES vehicles. As part of this effort, increasing TCA performance is one of 
the major design objectives for both engines. Original injector designs, dating to the 1960’s for the J-2 and 1990’s 
for the RS-68, were accomplished by the process noted above. Thus, the current designs are based on less than 
complete representation of the flow physics. Since essentially the same injector design tools are being used in the 
modification efforts today, the current development process relies on a similar mixture of empiricism and rules of 
thumb. This includes expensive subscale testing for design verification followed by the art of scaling these data to 
the full scale design. To effectively increase TCA performance beyond the current levels, improved design 
methodologies must be developed. 

Computational Fluid Dynamics (CFD), by replacing empiricism in the legacy design tools with a higher-fidelity 
physics-based approach, has the potential to improve the historical injector design and development process in at 
least three areas. First, initial designs can be improved by including relevant physical processes in the CFD models. 
Second, CFD models can be used to help create and execute an improved subscale test program. Third, the models 
can be used (whether in subscale or full scale testing) to interpret and understand test data. The physics-based 
approach is capable of including small injector geometry details that are known to have first order effects on key 
physical processes. It can also capture multi-dimensional and time dependent aspects of injector-driven flows, which 
are critical to achieving a robust design. Extensive parametric evaluations of the design space can be conducted to 
determine which variables are relevant and over what ranges. This information can in turn be used to decide the 
order and parametric combinations that should be tested to most effectively inform the full scale design. If there are 
schedule and budget constraints, the CFD model can also be used to augment test data so that a more complete 
understanding of detailed design issues can be achieved. An important use of a rigorously validated CFD model will 
be to provide unambiguous scaling of the subscale test data to the full scale design. 3 

Despite its potential, CFD is not currently employed in the injector design process to anywhere near the extent 
described above. This is largely due to the immaturity of the tool and significant uncertainties with respect to the 
models. From the perspective of injector design, CFD has not reached its potential for reasons that fall into three 
general categories. 4 The first is model fidelity, both in terms of geometry and physics. The second is solution 
robustness, i.e., the ability to achieve sufficiently accurate solutions in the timeframe dictated by the design process. 
The third is a lack of demonstrated accuracy of current solutions. Many times, especially in the production 
environment required for design support, injectors are modeled by assuming the flow is steady and axisymmetric. 
Similarly, propellants are often modeled as ideal gases when they are actually real fluids or in a liquid state. 
Continued increases in computer speed have elevated the accuracy of steady axisymmetric injector models to a 
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sufficient robustness level. However, as will be shown below, there are still significant accuracy problems with these 
solutions that decrease their utility. 

There has been very little methodical validation done for rocket engine injectors. This is explained in part by the 
fact that relevant data are difficult to acquire at the harsh operating conditions associated with rocket propulsion. 
Likewise, the requirements on CFD code capabilities and the required model suite to perform such calculations is 
equally daunting. The other reality is that meaningful validation, especially for rocket engine injector flows, requires 
a long-term commitment to tedious, methodical, and extremely time-consuming research. Without a concerted effort 
to demonstrate and improve injector model accuracy, the use of CFD for injector design will be infrequent and 
limited to qualitative trend analysis. 


II. Background 

The need to provide improved support to the TCA design process in general, and injector design in particular, 
has led to significant efforts in CFD code development, verification, and validation at NASA Marshall Space Fight 
Center (MSFC). 5 A notable achievement has been the evaluation of the current production capability based on 
steady axisymmetric Reynolds-Averaged Navier-Stokes (RANS) simulations. A primary application of this 
capability has been focused on shear coaxial elements at the single element model problem level. A series of 
experiments testing a range of coaxial injector elements, both shear and swirl, using hydrogen and oxygen have been 
carried out and results have been reported. 6 ' 8 These single element model problems have been used to evaluate CFD 
model accuracy by MSFC and others. The salient features of these efforts are provided here for context. 

A. Experimental Description 

Single element injectors were tested in a cylindrical copper heat sink chamber with the hot wall heavily 
instrumented with coaxial thermocouples and Gardon heat flux gauges. The instrumentation produced a set of 
temperatures and corresponding heat fluxes at discrete axial locations along the chamber wall. These wall 
temperatures were used as boundary conditions in the CFD model. Numerous injector designs, chamber pressures, 
and propellant combinations were tested. 9,10 

B. Demonstrated Accuracy Levels for Coaxial Single Element Injectors 

Ten of the experimental cases were selected to demonstrate CFD model accuracy. Each element type and 
propellant combination selected is representative of those used in actual rocket injectors. These three coaxial 
elements have served to establish the current RANS capability for treatment of single elements, which provides the 
starting point for the current effort. Element 1 operates with gaseous propellants flowing axially in concentric tubes 
(i.e., shear-coaxial injection). Element 2 is similar in geometry to Element 1 and retains the axial propellant flow. 
For this case, however, the central oxidizer jet is in a liquid state. Alternatively, Element 3 introduces the liquid 
oxidizer tangentially, imparting a swirling component to the flow. Thus, from Element 1 to Element 3 there is a 
progression of complexity in terms of geometry and flow physics. These three target cases were used to assess the 
capabilities of the current RANS CFD modeling approach. 

A series of calculations were performed to simulate the wall heat flux over a range of conditions and the results 
were compared to the experimental data. There were two specific objectives. The first was to evaluate the ability of 
the axisymmetric RANS CFD model to simulate the heat flux on the chamber wall produced by the reacting 
propellants. The second was to implicitly evaluate the model’s ability to predict the rate at which the propellants 
mixed and burned. This was to be inferred by the rise in the heat flux rate in the near injector portion of the chamber. 
A representative comparison of the results with the corresponding experimental data is shown in Figure 1 . 
Evaluating the solution accuracy was facilitated by dividing the flow domain axially into a near-injector region 
(where propellants mix and bum) and a downstream region (where burning is complete and combustion products 
mix with the remaining fuel). These regions are defined as the “head end,” which covers an interval 2 inches (5.08 
cm) from the injector and “downstream,” which is 8 inches (20.32 cm) from the injector, respectively. 

Accuracy levels were determined based on a zero to five range from the MSFC Simulation Readiness Level 
(SRL) credibility scale. Details of the SRL scale can be found in Reference 7. Evaluation of the accuracy level for 
each calculation is shown in Table 1. Examination of Figure 1 and Table 1 show that the data comparison is fair for 
Element 1 but is degraded as more physics and geometric complexity are added. The demonstrated accuracy level of 
all three calculations requires improvement before the RANS model can be used to quantitatively impact the injector 
design process. Achieving this goal is the motivation for the current effort. 

C. Paths to Improved Levels of Demonstrated Accuracy 
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One of the most significant shortcomings with attempts at validation is that there is no guidance on a path to 
improvements in accuracy. The current discrepancies between the RANS solutions and experimental data may be a 
product of faulty prediction of the propellant mixing and/or combustion characteristics. Errors could also lie in the 
model’s ability to predict heat fluxes, even with correct gas temperatures. Unfortunately, the precise cause (or 
causes) cannot be established definitively from the current set of validation cases and no definitive conclusions can 
be drawn regarding potential model improvements. Since conclusions are required to move the larger injector 
validation effort forward, the return on the resources invested previously has not been as high as desired. 

Philosophically, there are two approaches that can potentially lead to more concrete conclusions toward an 
advanced predictive capability. The first is based on a hierarchical approach . 11 The key element to this approach is 
the decomposition of a very large complex system, such as the TCA, into combinations of simpler problems. The 
injector, for example, is a TCA subsystem that can be decomposed into a set of model problems such as the uni- 
element rocket configuration . 4 5 The model problems are then further decomposed into unit physics problems, with 
each focused on isolated phenomena. Each unit physics problem employs a simple geometry that supports detailed 
measurements and uncertainty analysis. Using this building block approach, methodical validation of a particular 
CFD technique can potentially be obtained through definitive conclusions at each step. When there are no resource 
constraints, such a logical step-wise approach is the preferred validation path. In practice, however, execution of this 
approach is limited by constraints in both time and budget. The injector problem itself, for example, decomposes 
into a large set of unit physics problems. Applying an entire set of unit physics problems to validate a code is 
effectively impossible and thus requires prioritization. Consideration of only a subset of the unit physics problems 
limits coverage at this fundamental level and reduces the probability of success. 

A second approach to that described above is to use the existing single element heat flux data as a common 
reference and apply several different computational techniques with different levels of fidelity. Using this approach, 
the goal is to determine what level of accuracy can be achieved with a hierarchy of techniques applied to common 
“target cases” using identical geometric, boundary, and initial conditions. Higher fidelity methodologies, such as the 
Large Eddy Simulation (LES) technique, relax the RANS assumption of steady axisymmetric flow and presumably 
provide a higher level of accuracy. There are many issues, however, that must be addressed before such conclusions 
can be made with certainty. Modeling errors, numerical errors, and the validity of the fundamental assumptions that 
each technique is based on must be assessed relative to respective target cases before quantitative advances in 
predictive capabilities can be made. In particular, all simulations incorporate some form of formal filtering of the 
governing conservation equations. The RANS and unsteady-RANS (URANS) approximations have historically 
provided the central basis for engineering design codes throughout industry and government. RANS employs 
filtering in time to derive the governing conservation equations for the mean state. All dynamic degrees of freedom 
smaller than the largest energy containing eddies are averaged and no information exists to describe unsteady 
interactions between the small-scales. The primary advantage of RANS is that it requires a minimum of 
computational resources, which facilitates the fast solution turn-around times required for most engineering design 
studies. A disadvantage of RANS, however, is that all of the turbulent motions are modeled and thus it only provides 
a bulk representation of the strongly coupled velocity, scalar mixing and chemical interactions that occur over an 
extremely wide range of turbulence scales. 

In contrast to RANS, LES has historically employed spatial filtering to split the field variables into time- 
dependent resolved-scale and subgrid-scale (SGS) components. The large energetic scales are resolved and the SGS 
quantities are modeled to provide a complete time-accurate representation of dynamic processes over the full range 
of relevant time and length scales. The potential advantage of LES over RANS (and URANS) is due to the fact that, 
when properly implemented, it provides a complete multi-scale closure that describes the full range of dynamically 
active scales in a given system. RANS by design, on the other hand, only provides a bulk representation of these 
dynamic interactions. The disadvantage of LES over RANS, of course, is the cost. Algorithmic requirements must 
also be considered. Compared to RANS, the use of LES imposes an extremely strict set of algorithmic and 
resolution requirements that must be rigorously enforced. A significant range of energetic scales must be resolved 
with minimal computational errors. This has stringent implications regarding grid quality and the amount of 
allowable grid stretching that can be used. Numerical dissipation and dispersion errors must also be minimized since 
they have significantly detrimental effects on SGS models. The presence of these errors depletes energetic 
turbulence scales at the mid- to high-wavenumbers and consequently competes with the models. When this occurs 
the SGS models themselves have little or no effect and the contamination often leads to erroneous conclusions. 

Regardless of the validation approach and CFD technique used, understanding the tradeoffs listed above are 
central to the development of an improved level of demonstrated accuracy. Each of the validation approaches has 
advantages and disadvantages and neither is guaranteed to be successful. The first is based on a hierarchical 
approach where complex systems, such as the TCA, are decomposed into combinations of simpler unit problems. 
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The second applies a hierarchy of CFD techniques to common TCA relevant target cases using identical geometric, 
boundary, and initial conditions with the goal of determining what level of accuracy can be achieved, and more 
importantly, how higher fidelity (albeit more expensive) methods can be used to improve the accuracy of lower 
fidelity (albeit faster and more economical) engineering approaches. Considering the need to support Constellation 
propulsion design projects with more quantitative analyses, the decision has been made to proceed along both paths 
simultaneously. 


III. Scope of the Current Effort 

The current validation effort applies the second of the two approaches described above. We have adopted a 
common set of target cases, which includes; the shear coaxial element with gaseous oxygen (G0 2 ) and gaseous 
hydrogen (GH 2 ) propellants, the shear coaxial element with liquid oxygen (L0 2 ) and GH 2 propellants, and the swirl 
coaxial element with G0 2 and GH 2 propellants. Each of these three target cases will be simulated with a group of 
tools that covers a range of simulation fidelities, including RANS, URANS, and LES. The goal is to determine what 
degree of demonstrated accuracy is possible at each level of fidelity and the cost, in terms of time and manpower, 
required to achieve these different levels of accuracy. Put in the context of injector design support, the goal is to 
understand what levels of demonstrated accuracy are possible and what levels are affordable in terms of resources in 
the production environment. This paper documents our current status with respect to the first model problem; the 
single element shear coaxial injector with gaseous oxygen and gaseous hydrogen propellants. Our current focus is 
the initial evaluation arid comparison of the results in the near injector flow field where the gases begin to mix and 
react. This region is largely responsible for the flame dynamics that govern the rest of the flow field. Comparison of 
the results to the experimental heat flux data will be presented in a future paper. 

A. Description of the Experiment 

The experiment that provided the data for Element 1 was conducted in the Cryogenic Combustion Laboratory at 
The Pennsylvania State University. The test rig, shown schematically in Figure 2, is comprised of oxidizer and fuel 
prebumers, a single element shear coaxial injector, and an instrumented heat sink main combustion chamber. The 
experiment was designed to measure the axial heat flux profile along the chamber wall and a corresponding wall 
temperature profile to be used as a boundary condition for the simulations. Coaxial heat flux gauges in the chamber 
wall provide data for both the wall temperature and corresponding wall heat flux profiles. The test simulated in this 
effort was conducted for a target chamber pressure of 5.5MPa. The propellant mixture ratio was 6.7. The test rig, 
instrumentation, flow conditions and data are described in detail in Reference 9. 

The computational domain and common set of boundary conditions used in the calculations are shown in Figure 
3. The overall domain is cylindrical, has been non-dimensionalized using the oxidizer post inner diameter (5 = 5.26 
mm) as the reference length scale, and includes the injector, main chamber and nozzle. All three sections match the 
geometric profiles of the actual experimental apparatus. The injector section is 28.973 units long, which provides 
appropriate entry lengths required for development of the turbulent boundary layers in the fuel and oxidizer streams. 
The total length of the combustion chamber and nozzle is 63.983 units (54.325 and 9.658 units, respectively). The 
inlet diameter of the central oxidizer jet is 1.5057 units. The exit diameter is 1 unit, and the oxidizer post is recessed 
0.081749 units behind the chamber face. The inner and outer diameters at the inlet of the annular fuel jet are 2.4144 
and 4.8288 units, respectively. The corresponding exit diameters are 1.1977 and 1.4240, respectively. The diameters 
of the combustion chamber, nozzle throat, and nozzle exit are 7.2434, 1.5525, and 2.2750, respectively. 

The composition, key reference properties, and flow characteristics of the oxidizer and fuel streams are listed in 
Table 2. The oxidizer stream contains 0.906 moles of oxygen, 0.0940 moles of water, and is injected at a 
temperature of 711 K. The fuel stream contains 0.857 moles of hydrogen, 0.143 moles of water, and is injected at 
800 K. Calculations were performed by applying constant flow rates of 0.0904 kg/s (oxidizer) and 0.033 1 kg/s (fuel) 
at respective inlets. Supersonic boundary conditions are applied at the nozzle exit. No-slip conditions are applied at 
all wall surfaces. All of the injector wall surfaces are assumed to be adiabatic. The face of the combustion chamber 
and oxidizer post tip are assumed to be fixed at 755 K. The axial temperature distribution along the radial wall of the 
combustion chamber is obtained using available experimental measurements. A least squares fit of these data is used 
to approximate this distribution, as shown in Figure 2. The nozzle temperature is assumed to be fixed at 510 K. This 
value is consistent with the curve-fit, which by design provides a value of 510 K at the end of the combustion 
chamber with a gradient of zero. The reference pressure for the current calculations was set at 5.2 MPa. At this 
pressure the bulk exit velocities of the injector are 154 m/s for the oxidizer stream and 764 m/s for the fuel stream. 
The corresponding Reynolds numbers based on the hydraulic diameter (i.e., 8 for the oxidizer stream and d 0 - dj for 
the fuel stream) are 604,000 and 169,000, respectively. 
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B. General Approach 

Using the common configuration and conditions described in Figure 3 and Table 2, we have applied a 
progressive hierarchy of 5 computational techniques. Here and in the results we will refer to the respective 
approaches using the co-authors names - Oefelein: LES using a Stochastic Reconstruction Model, Menon: LES 
using the Linear Eddy Model, Yang: LES with a Laminar Flamelet Model, Merkle: Unsteady RANS, and Tucker: 
Steady RANS. Note that Oefelein and Menon have performed three-dimensional calculations. Yang, Merkle and 
Tucker have performed two-dimensional axisymmetric calculations. The level of spatial/temporal fidelity is 
systematically increased as we move from Steady RANS to high-fidelity LES. Naturally, the cost also increases 
proportionately. Subsections C-G provide a description of each approach along with pertinent details related to the 
implementation. 

C. LES - Stochastic Reconstruction Model (Oefelein, Sandia) 

Computational Model: The baseline theoretical-numerical framework used for this method solves the fully 
coupled conservation equations of mass, momentum, total-energy, and species for a general class of chemically 
reacting flows. The overall model handles both multicomponent and mixture-averaged systems, with a generalized 
treatment of the equation of state, thermodynamics, and transport processes. Details are given in References 12-14. 
Here we employ a new class of reconstruction models 15,16 that combines the purely mathematical approximate 
deconvolution procedure with physical information from an assumed scalar spectrum to match specific scalar 
moments. Using this method, a surrogate to the exact scalar field can be estimated such that the filtered moments 
match to a specified order. In principle, the surrogate field can be used to calculate the SGS contribution of any 
related nonlinear function. In practice, however, the extent of the nonlinearity limits the accuracy and it has been 
shown that this method cannot be used reliably to close the filtered chemical source terms directly. It can be used, on 
the other hand, to obtain highly accurate representations of polynomial noniinearities such as the SGS scalar 
variances. These are precisely the input required to generate SGS fluctuations stochastically. Given these findings, 
we have proposed an extension to this approach by coupling it to a stochastic reconstruction model. The combined 
methodology provides a correlated approximation of the SGS velocity and scalar fluctuations with the correct time 
history and spatial distribution. The modeled instantaneous field (i.e., (p = <<p> + <p', where «p> represents the 
resolved-scale contribution of an arbitrary scalar and cp' the correlated SGS fluctuation) is used to evaluate the 
filtered finite-rate chemical source terms directly. The filtered source terms, and related scalar field, are closed by 
selecting an appropriate detailed chemical kinetics mechanism in the same manner as is done for a Direct Numerical 
Simulation (DNS). The model coefficients are evaluated locally in a manner consistent with the dynamic modeling 
procedure. Thus, the only adjustable parameters are the grid spacing, integration time-step, and boundary conditions. 
In the limit as the grid resolution and time-step approach the smallest relevant scales, contributions from the SGS 
models (i.e., the fluctuating component) approach zero and the solution becomes a DNS. 

Model Implementation and Baseline Grid : The numerical framework has been optimized to meet the strict 
algorithmic requirements imposed by LES. The governing conservation equations are integrated in time using the 
dual-time stepping technique, with a generalized preconditioning methodology that treats convective, diffusive, 
geometric, and source term anomalies in an optimal manner. The spatial scheme employs staggered finite-volume 
differencing in generalized curvilinear coordinates. This approach provides non-dissipative spectrally clean damping 
characteristics and discrete conservation of mass, momentum and total-energy, which is a critically important 
feature for LES. The algorithm has been optimized to provide excellent parallel scalability attributes using 
distributed multiblock domain decomposition with generalized connectivity. Distributed-memory message-passing 
is performed using MPI. The software package has been ported to all contemporary massively-parallel computer 
platforms and is well suited for high-performance computations. The calculation has been performed in two stages. 
First, the injector section was run to establish the time-dependent flow condition at x/8 = -5. This was then used as 
the inflow condition for the main calculation, which is currently being run on 528 processors with full chemistry, 
detailed thermodynamics and full transport. Finite-rate chemical kinetics are modeled using the hydrogen-oxygen 
mechanism developed and optimized by Conaire et al. 17 The baseline grid resolution is established at the injector 
exit (x/5 = 0, see Figure 3) by insuring the grid spacing in wall units is consistent in all three directions and small 
enough to resolve the boundary layer profile and stresses (i.e., the wall-resolved approach). Resolving this level of 
energy has a direct impact on the accuracy with which the shear-layers downstream of the injector evolve. The most 
desirable spacing is obtained by stretching in the wall-normal direction such that the first cell from the wall is within 
a y + value of 1, with the first 16 cells within the interval 0 < y + <30 and Ay + « 50 in the core region. Similarly, the 
most desirable transverse grid spacing is Ax + « Az + ~ 50. Given the relatively high Reynolds numbers, applying the 
criteria outlined above in the strict sense is prohibitive. Thus, the grid has been optimized by maximizing the relative 
spacing in key sections while insuring that the near-wall boundary layer dynamics and peaks in the Reynolds-stress 
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tensor are being resolved to an acceptable level. The initial baseline grid used for this study is approximately 33- 
million cells. The (normalized) injector velocity profiles obtained using this level of resolution is shown in Figure 4. 
Contours represent the axial velocity component, vectors represent the transverse components. Note that the bulk 
velocity of the fuel stream is 4.96 times greater than the oxidizer stream. 

D. LES - Linear Eddy Model (Menon, G. Tech) 

Computational Model: For this case, the SGS closure employs a transport model for the sub-grid kinetic energy 
and a SGS EBU closure for the combustion process. 18,19 A 21-step; 8-species mechanism is employed assuming a 
thermally perfect gas. Simulation of the inlet is limited to a region 5 cm upstream of the dump plane and the bulk 
inflow conditions were chosen based on an independent simulation of the full injector region. Characteristic inflow 
and supersonic outflow are employed. Isothermal combustor walls based on experimental data and isothermal 
conditions are used along the injector walls. On our Intel Xeon cluster, it takes about 5200 single-processor hours to 
simulate about 0.5 ms of physical time. Approximately 8 ms of data has been time-averaged after the initial 
transients. 

Model Implementation and Baseline Grid : The simulation software LESLIE3D is a finite-volume solver 
developed at Georgia Tech that is second order accurate in space and time. It employs a hybrid approach that 
combines a predictor-corrector scheme in smooth flow regions and a MUSCL flux extrapolation with a Harten-Lax- 
Leer (HLL) type limiter in regions of very high gradients and/or contact discontinuities. No explicit dissipation is 
employed. Switching between the two algorithms is dynamic and local. 20 The grid is cylindrical with dimensions of 
611x87x65 and an inner “butterfly” Cartesian section of dimensions 611x17x17 to eliminate the singularity. 
Clustering is employed in the shear layer and in the post region. 

E. LES - Flamelet Model (Yang, PSU) 

Computational Model: The basis of this approach is the general theoretical/numerical framework described in 
References 21-23. The formulation accommodates the full conservation laws for a multicomponent chemically 
reacting system. Full account is taken of real-fluid thermodynamics and transport over the entire temperature and 
pressure regimes of concern. Thermodynamic properties, such as internal energy, enthalpy, and constant-pressure 
specific heat, are obtained directly from fundamental thermodynamics theories. 24 Transport properties are estimated 
using an extended corresponding-state theory. SGS turbulence-chemistry interactions are treated by means of a 
mixture-fraction based flamelet model. The approach assumes that chemical scales are shorter than the Kolmogorov 
scale of turbulent flows. Consequently, a turbulent flame can be envisioned as a synthesis of thin reaction zones 
embedded in an otherwise inert turbulent flowfield and the inner structure of the flame can be handled separately 
from turbulent flow simulations. Instead of directly treating reactive scalars, the focus of the flamelet model is 
placed on the identification of the flame surface in the flowfield, which can be obtained by solving the conservation 
equation of the mixture fraction in a coupled manner with the mass, momentum, and energy equations. Because the 
flame thickness is typically smaller than the grid size, the influence of the SGS mixture fraction variance on the 
flame behavior is modeled either by a presumed 5- or a P-shaped probability density function (PDF). Once the 
solution of the mixture fraction and its variance are acquired, the species mass fractions are evaluated using a pre- 
established flamelet library based on a series of numerical solutions of counterflow hydrogen/oxygen diffusion 
flames. 25 The flame solutions were carried over a wide range of flow strain rates with a detailed reaction mechanism 
including 8 species (i.e., H 2 , 0 2 , H, O, OH, H0 2 , H 2 0, and H 2 0 2 ) and 19 reversible reactions. The validity of the 
flamelet model under the present simulation conditions has been confirmed through a careful comparison of the 
characteristic chemical and turbulent time scales over the entire flowfield. 26 Specifically, the Kolmogorov time scale 
evaluated based on the LES simulation results is 1.2 x 10" 6 s at the axial location 5 mm downstream of the injector 
faceplate. The characteristic chemical time scale of 5.48 x 10' 7 s obtained from flame calculation is at least 2 times 
smaller than this value. 

Model Implementation and Baseline Grid : The simulation was initialized using the RANS approximation to 
establish a reasonable initial condition for the LES calculation. The combustor was preconditioned with a mixture of 
oxygen (0.3 by mass), hydrogen (0.3), and water vapor (0.4). The chamber pressure was set at 5.2 MPa. The steady- 
state solution obtained from the RANS calculation was then interpolated and mapped onto the LES grid as the initial 
condition for the LES simulation. The numerical implementation includes the incorporation of general-fluid 
thermodynamics and transport theories into a preconditioning scheme along with a dual time-stepping integration 
algorithm. 27,28 All the numerical properties, including the preconditioning matrix, Jacobian matrices, and 
eigenvalues, are derived directly from fundamental thermodynamics theories, rendering a self-consistent and robust 
algorithm. A multiblock domain decomposition technique, along with static load balance, is employed to facilitate 
the implementation of a distributed parallel computation using MPI. The parallelization methodology is robust and 
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the speedup is almost linear. A grid of 750x320 cells in the axial and radial directions, respectively, was employed. 
Special attention was paid to the grid resolution close to the chamber wall in order to ensure an accurate prediction 
of wall heat flux. The grid point next to the wall is 10 pm. The initial grid stretching ratio is 1.06. The computation 
domain was divided into 71 blocks, with each calculated on a 2.2 GHz Pentium IV processor of a distributed- 
memory parallel computer. The physical time step was 1.0 x 10" 7 s, and the maximum CFL number for the inner 
loop integration in pseudo-time was 0.7. The simulation was conducted for an extended period until the flowfield 
reached its stationary state. Data was then collected for 4 flow-through times (i.e., 5 ms) to obtain statistically 
meaningful mean properties. The overall computation took approximately 60,000 CPU hours. 

F. Unsteady RANS (Merkle, Purdue) 

Computational Model: The Spalart-Allmaras turbulence model was used with a 9 species, 17 reaction step model 
for hydrogen and oxygen. Enthalpies, viscosity and thermal conductivity for each species are taken from polynomial 
expressions. Species diffusivities are obtained from a Chapman-Enskog expression. Additional details are given in 
References 29 and 30. The computations are based on a two-dimensional (axisymmetric) approximation with flow 
conditions analogous to those presented in Table 2. 

Model Implementation and Baseline Grid : The URANS code employed in this study uses an implicit, second- 
order accurate scheme based upon dual time stepping and a modified Roe-based flux scheme that controls artificial 
dissipation in transients and low-speed regions. A total of approximately 200,000 cells were used for the 
simulations. As initial conditions for the computation, the first three-fourths of the inlet channels were filled with 
fuel and oxidizer respectively at the appropriate temperatures given in Table 2. The remainder of the two inlet 
channels and the entire combustor were filled with nitrogen at 1500 K. Details are provided in Section IV below and 
contrasted with the other approaches. The initial nitrogen in the chamber and inlet tubes served two purposes. First, 
since neither the incoming fuel nor oxygen contain nitrogen, the amount of residual N 2 in the computational domain 
immediately indicates the degree to which the initial condition still impacts the predictions. Unambiguous results 
clearly require that the nitrogen concentration be reduced to a negligible magnitude at every point throughout the 
domain. The second function of the nitrogen was to prevent combustion between hydrogen and oxygen until mixing 
started, with the 1500AT initial temperature being chosen to ensure spontaneous ignition as soon as fuel and oxidizer 
came together so that no amount of pre-mixing would be tolerated. This resulted in a very smooth and well 
controlled ignition process that avoided large acoustic oscillations caused by ignition of a finite mass of propellant. 

G. Steady RANS (Tucker, MSFC) 

Computational Model: The Mentor Shear Stress Transport 31 (SST) and original Wilcox 32 k-omega models were 
used in the current study. The model used for finite-rate hydrogen-oxygen chemistry was the 6 species 28 reaction 
scheme. 33 Thermodynamic properties are obtained using a standard partition function formulation which calculates 
the specific heats, internal energies and entropies of each individual perfect gas species. A typical solution requires 
about 10,000 time steps to achieve steady state convergence. A typical calculation using 16 AMD Opteron 246 
processors requires approximately 4 days of wall clock time or 1,600 CPU hours. As long as the number of cells per 
processor remains similar to that described above, the Loci-CHEM CFD program (as described below) is nearly 
perfectly scalable. 

Model Implementation and Baseline Grid : Loci-CHEM is a finite-volume flow solver for generalized grids 
developed at Mississippi State University in part through NASA and NSF funded efforts. CHEM uses high 
resolution approximate Riemann solvers to solve turbulent flows with finite-rate chemistry. Preconditioning 34 is 
available for low Mach number applications. Details of the numerical formulation are given in the CHEM user 
guide. 35 Loci-CHEM is comprised entirely of C and C++ code and is supported on all popular UNIX variants and 
compilers. Parallelism is supplied by the Loci 36 framework which exploits multi-threaded and MPI libraries. All 
meshes are generated using GRIDGEN Version 15 from Pointwise Corporatipn. 37 Hybrid meshes were generated by 
extruding connecters representing the axisymmetric solid surfaces normal to the solid surface into the computational 
domain with an initial spacing and stretch rate defined. The extrusion was terminated at the location where the 
extruded cells possessed an aspect ratio of approximately unity. This allowed a good quality transition from 
structured to unstructured cell types. The remainder of the computational domain was meshed using unstructured 
cells with a boundary decay factor of approximately 0.995. Inlet boundary values for turbulence quantities, k and 
omega, were specified as 5xl0" 6 m 2 /s 2 and 500 1/s respectively. While these values are lower than those 
recommended, 31 the turbulent field develops in the injector inlets and previous experience has shown that the 
injector inlet distances in this case allow no sensitivity of the chamber solution to the injector inlet turbulence 
quantities specification. The computational domain was initialized with quiescent steam at 780 K to start the steady 
state simulation. The steam provided a high enough initial temperature such that when the reactant streams reached 
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the post-tip region and mixed self-ignition occurred. The model was implemented in local time stepping mode with 
a physical time step of lxl O' 4 seconds, which was dynamically reduced locally to not exceed a maximum CFL 
number of 10,000 or a change in temperature, pressure or density variables of more than 10 percent. With these 
settings, the cases ignited relatively smoothly and ignition did not cause any sustained reverse flow of propellants 
into either the fuel or oxidizer inlet tubes. Solution convergence was evaluated by a combination of residual drop, 
total mass flow convergence in terms of integrated mass at inlet verses integrated mass at outlet, species stored mass 
convergence integrated over the entire domain volume, and temperature and pressure iteration history at selected 
probe point locations throughout the flow-field. In addition the chamber wall heat flux was monitored until 
convergence was identified. Typically three orders of magnitude of residual drop are obtained. Computational 
domain mass conservation was always achieved to with 0.1 percent of the total mass flow. The species stored 
masses are also driven down to less than 0.1 percent of the total stored mass of that species. Similarly, the probe 
point pressures and temperatures are also driven down to less than 0.1 percent. 

IV. Results 

Our current objectives are to establish the initial baseline observations, to sort out and highlight the critical 
differences in the way the calculations are being done, to comment on the potential impact on the interpretation of 
results, then to suggest improved "best practices" for future calculations. Where prudent, we provide a list of 
observations regarding the differences and similarities in each solution with potential pros/cons affecting solution 
accuracy. We have stopped short on any speculation at this stage and primarily focus on clearly outlining the 
observations, with suggestions on how to improve the quality of future calculations. 

An issue of current interest is the way in which each of the calculations has been initialized and the potential 
impact on the observed results. Of the five methodologies described in Section III, the calculations performed by 
Oefelein, Menon, Yang and Merkle are unsteady calculations. The calculation performed by Tucker represents the 
baseline steady RANS calculation. Three-dimensional (3D) calculations were performed by Oefelein and Menon. 
The remaining calculations are two-dimensional (2D) axisymmetric. From the fluid dynamic perspective, each of 
the different modeling approaches was applied using a common set of boundary conditions and there was good 
consistency between the cases. There are many interesting differences, however, in the way this initial set of 
calculations was initialized and a wide variety of grid resolutions were used. Different chemical kinetics models are 
also being used at present. 

By definition, each modeling approach (and associated numerical scheme) has a different set of spatial and 
temporal resolution requirements. The initialization procedure, however, is something that should be implemented in 
a fairly uniform manner, independent of the method used, and should also have a well defined criteria associated 
with it. In the discussion that follows, we define a flow through time as that required for a particle to traverse the 
entire length of the chamber (i.e., from the injector exit plane to the nozzle exit) at a bulk velocity based on the total 
mass flow entering the system and the cross-sectional area. Assuming theoretical products at 3410 K yields a bulk 
chamber velocity of approximately 41 m/s, which translates to a flow through time of 8.3 ms. Using this definition, 
initialization of respective fields downstream of the injector for the current set of results was performed as follows: 

• Oefelein: Theoretical product mixture using a flamelet solution at a mixture fraction of 0.268. The volumetric 
composition and temperature of the product field were 18.6% (H 2 ), 0.74% (0 2 ), 5.58% (OH), 75.08% (H 2 0), 
and 3410 K. The initial reference pressure was set to 5.2 MPa. The solution was initialized by running the full 
case on a 2D radial sector of the 3D grid for 4 flow through times, then mapping the solution to the full 3D 
grid and running for an additional 2 flow through times. 

• Menon: Approximated the product mixture using a volumetric composition and temperature of 90% (H 2 0), 
10% (OH), and 3000 K. The initial reference pressure was set to 5.2 MPa. The calculation was initially started 
with a reduced 2-step mechanism and ran for approximately 15 ms to flush out the initial transient. After this 
the full mechanism was employed and ran about 5 ms. 

• Yang: Initialized the calculations using a steady RANS solution to establish a reasonable initial condition for 
the LES calculation. The combustor was preconditioned with a mixture of 30% (H 2 ), 30% (0 2 ) and 40% H 2 0. 
The chamber pressure was set at 5.2 MPa. The steady-state solution obtained from the RANS calculation was 
then interpolated and mapped onto the LES grid. The LES calculation then ran for 4 ms to allow the flow field 
to reach its stationary state. 

• Merkle: Applied an initialization procedure where the first three-fourths of the inlet channels were filled with 
fuel and oxidizer respectively at the appropriate temperatures. The remainder of the two inlet channels and the 
entire combustor were filled with nitrogen at 1500 K. The initial pressure was set to 3.24 MPa, which 
corresponds to pressure generated in the chamber with 1500 K nitrogen passing through the choked throat at 
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the same mass flow rate as the incoming propellants. As the fluid exiting the choked throat increased in 
temperature, the chamber pressure gradually increased to near the observed experimental condition. The initial 
velocity throughout the entire chamber was set to zero with an ensuing shock used to initiate flow. 

• Tucker: The steady RANS calculation was initialized with quiescent steam at 780 K to start the steady state 
simulation. The steam provided a high enough initial temperature such that when the reactant streams reached 
the post-tip region and mixed self-ignition occurred. The model was implemented in local time stepping mode 
with a physical time step of lxl 0" 4 seconds, which was dynamically reduced locally to not exceed a maximum 
CFL number of 10,000 or a change in temperature, pressure or density variables of more than 10 percent. With 
these settings, the cases ignited relatively smoothly and ignition did not cause any sustained reverse flow of 
propellants into either the fuel or oxidizer inlet tubes. 

To investigate the preliminary differences between these solutions, we first present a set of key instantaneous 
results, then the corresponding time-averaged (and steady-state) solutions. It is important to point out that this is the 
first initial set of results. The primary goal here is to better understand how to unify our basis for comparison, to 
improve the consistency with which the calculations are performed, and to establish a consistent criterion for future 
calculations. Our initial focal point is a radial cross-section in the near injector region which spans the ranges of 
-2 < x < 22 mm and 0 < r < 6 mm. 

General Flowfield Characteristics: 

Figure 5a-e show the instantaneous fields of temperature, vorticity magnitude, and mass fractions of hydrogen, 
oxygen and OH radical in the vicinity of the oxidizer post for each of the 4 unsteady calculations. As is typical for 
this configuration, a diffusion-dominated flame emanates immediately downstream from the post and propagates 
along the boundary of oxidizer stream. A wake region, which consists of hot combustion products, effectively 
separates the fuel and oxidizer streams. The near-field flow dynamics are characterized by the evolution of the three 
mixing layers originating from the inner and outer edges of the fuel annulus and the inner rim of the oxidizer post. 
The development of the inner mixing layer surrounding the oxidizer, however, is inhibited by the expansion of the 
combustion products in the flame zone. Energetic vortices emerging from the oxidizer post evolve in a complex 
manner within the small recirculation zone just downstream. 

Figure 6a-f shows the corresponding time-averaged distributions. Here we compare temperature, vorticity 
magnitude, turbulent kinetic energy (TKE), and mass fractions of hydrogen, oxygen and OH radical in the near- 
field. There are typically two high temperature regions observed in the mean field. A small intensive reaction zone 
appears immediately behind the oxidizer post and the main flame starts approximately 6 mm downstream. The 
existence of a relatively low-temperature region between these reaction zones may be attributed to the evolution of 
large-scale vortices in the near-field and the associated high flow strain-rates. The magnitude of the vorticity in the 
shear layers surrounding the fuel stream is typically much greater than that in the rest of the domain. Two intense 
TKE regions appear in the upper and lower edges of the fuel annulus and merge together approximately 2 to 3 mm 
downstream of the injector faceplate. The potential cores of the fuel and oxidizer streams are clearly observed in the 
hydrogen and oxygen mass fraction fields, respectively. 

Preliminary Observations: 

To compare the instantaneous characteristics of respective solutions, we consider the unsteady results from 
Oefelein, Menon, Yang and Merkle. Figure 5a shows the instantaneous temperature field from the four different 
solutions. Note that respective plots have been generated using the same color maps and contour ranges. Several 
interesting observations can be made. First, there is a somewhat wide disparity between the results. Yang and 
Menon, for example, show similar temperatures in the upper-left comer where the flow recirculates from the upper 
radial portion of the combustion chamber. Oefelein’s and Merkle’s solutions, on the other hand, are both higher in 
this region and comparable in magnitude. Interestingly, Oefelein started with a relatively high initial temperature in 
the chamber (i.e., that associated with the theoretical products) and Merkle started with the coldest initial 
temperature. This points directly to the initialization procedure and the total time the cases were run to eliminate any 
residual initial conditions as potential issues. There is also a distinct difference between the 2D structures (Yang and 
Merkle) and 3D structures (Oefelein and Menon). The 2D structures look qualitatively similar. The 3D structures, on 
the other hand, are somewhat different. 

Comparing the vorticity magnitude (Figure 5b) provides additional insights. Trends are similar to those observed 
in Figure 5a. Oefelein’s solution has a much finer vorticity distribution but is qualitatively similar to Yang’s and 
Merkle’s. Menon’s solution is more dispersed. The distributions of vorticity in the shear layer regions are also 
different between the cases. Vorticity from the annular hydrogen stream is much finer in Oefelein’s and Merkle’s 
solutions. They are somewhat thicker in Menon’s and Yang’s. Oefelein shows more detailed structures, as indicated 
by small vortices at the wall. 
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Vortical (i.e., convective) mixing clearly has a direct impact on respective mass fraction distributions. Note, for 
example, the hydrogen concentration in the upper-left comer of Figure 5c. There is a wide disparity in the 
concentration between the cases and it is not clear if this is a residual effect from the initialization procedure or due 
to something else. It is also interesting to note the difference in the structures. Comparing Oefelein and Merkle, the 
distributions are similar, but the structure is different, Merkle’s being more coherent. This may be attributed to the 
difference between the ring structure in 2D and the vortex stretching mechanism present in 3D calculations. There is 
also a distinct difference between the two 3D cases (Oefelein and Menon), possibly due to the differences in local 
grid resolution. 

The oxygen and OH radical distributions exhibit similar trends (Figures 5d and e). In Menon’s solution, both the 
oxygen and OH are more diffused than the others. The oxygen core is still intact for all the cases. Interestingly, the 
amount of oxygen in the small recirculation zone just downstream of the oxidizer post increases on going from 
Oefelein’s (top) to Merkle ’s (bottom) solution. This could be due to the difference in chemical mechanisms, among 
many other possibilities. In analyzing the qualitative trends associated with the OH radicals, it must be pointed out 
that a different chemical mechanism was used for each case. Oefelein’s solution has relatively higher concentrations 
of OH and a very fine structure. Differences in the magnitude of OH are striking and need to be investigated further. 

In Figure 6 we show a comparison between the time-averaged results. Here we compare solutions from Menon, 
Merkle, Yang and Tucker. Recall the first three are unsteady calculations, and Tucker’s solution is from a steady 
RANS calculation. Note that the unsteady cases were averaged in slightly different ways. Merkle, for example, 
sampled every 100 time steps rather than every step, so the residual structure is still visible. Menon and Yang, on the 
other hand, sampled much more frequently. Comparing the time-averaged solutions to the RANS-averaged solution 
shows many interesting differences. 

Figure 6a shows a comparison between the temperature fields. Menon’s and Merkle’s solutions are somewhat 
similar. Yang’s solution is closer to Tucker’s steady RANS. Interestingly, this is not precisely the same trend as was 
observed for the instantaneous solutions. All of the time averaged plots show an initial hot spot, a cooler region 
slightly downstream, and then a hot region again. The flame is most likely disrupted in the cooler region due to the 
high flow strain rates that exist there. The potential core associated with Menon’s and Yang’s solutions are relatively 
short. In terms of the overall structure, Menon’s and Merkle’s are the most similar. Also note that all three of the 
unsteady solutions are distinctly different from the steady RANS solution. Tucker’s solution also exhibits a much 
higher temperature than the three unsteady cases. This is most likely due to dynamic mixing processes that the 
unsteady calculations account for. By definition, there is no mechanism to provide this mixing in the steady RANS 
calculations. Thus, a key mixing mechanism is potentially missing that could have a profound effect on the heat 
transfer characteristics (among other things). This observation is corroborated by the trends shown in Figures 6b and 
c. The steady-RANS results (Tucker) exhibit a much lower level of vorticity relative to the unsteady calculations. 
This is indicative of a lower level of convective mixing. The turbulent kinetic energy distributions shown in Figure 
6c are similar, which suggests that there is also a lack of diffusive mixing. Interestingly, the three unsteady 
calculations are similar in the overall vorticity and TKE distributions. 

Analysis of the time-averaged mass fraction distributions suggests that the initial condition is extremely 
dominant. Distributions of hydrogen are shown in Figure 6d. There is a distinct disparity between each case. 
Menon’s and Merkle’s solutions are qualitatively similar. Yang’s solution exhibits a higher level of hydrogen in the 
upper-left comer. Tucker’s solution exhibits levels in between. Yang’s is above Tucker and potentially high due to 
residuals from the initial solution. Menon’s is lowest. Interestingly, Merkle’s 2D solution, which was initially started 
with no OH in the chamber, has higher relative concentrations than Menon’s 3D solution. The cause is not clear. 
Similar trends are observed in Figures 6e and f. Menon and Merkle’s solutions are similar. Yang and Tucker’s are 
somewhat similar except for the fact that the steady-RANS has a longer potential core. The steady RANS solution 
has a noticeably longer penetration length compared to all of the unsteady calculations. Oxygen is clearly consumed 
faster in the unsteady calculations, which is again indicative that the mixing mechanism(s) inherently missing in the 
steady RANS calculations are important. 


V, Conclusions and Future Work 

The results presented above have provided a first glimpse into the many intricacies associated with establishing a 
predictive capability for rocket injector design. The largest uncertainty with the current calculations is the manner in 
which the solutions are being initialized and the amount of time required to eliminate residuals from the initial 
condition. Similar chemical mechanisms must also be used to eliminate ambiguities associated with species 
distributions and. the overall heat release rates. For the latter, we have proposed that the Conair mechanism 
(Reference 1 7) be used since this is the most recent. A clear set of spatial and temporal resolution requirements must 
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be established. It will also be useful to look at both smaller domains, such as the wake region immediately 
downstream of the oxidizer post, as well as the entire chamber to better understand how small-scale dynamics in 
critical regions of the flow affect the larger-scale processes and overall system performance. 

After establishing a clear criterion for implementation, systematic comparisons between the same hierarchies of 
modeling techniques must be made that include analysis of key instantaneous and statistical quantities. This includes 
1) snapshots of the flow field and flame structures in critical regions of the flow, 2) spectral information at selected 
locations (i.e., shear layers and flame fronts) to quantify the flow and flame dynamics, and 3) time-averaged flow 
properties (i.e., mean and rms velocities, turbulence properties, temperature, and species concentrations). It will also 
become more and more important to analyze animations of key dynamic processes to better understand how to 
interpret various results and averaged data. Such analyses will serve as a precursor for making accurate predictions 
phenomena such as wall heat fluxes, which is one of the central objectives. 

As part of establishing a predictive capability, it will become increasingly important to prescribe accurate 
implementation procedures and equally important to understand the investment in resources required to achieve 
certain levels of accuracy. Future efforts will be aimed toward these objectives. 


Acknowledgments 

Co-workers at MSFC who executed the steady-RANS calculations are Jeff Lin, Jeff West and Doug Westra. 
Graduate Students who helped, etc., from G. Tech., Nan Zong (Research Associated at PSU), Purdue. NASA CUIP 
Program. NASA funding for JCO. 

References 

'Muss, J.A., Nguyen, T.V., and Johnson, C.W., "User’s Manual for Rocket Combustor Interactive Design (ROCCID) and 
Analysis Computer Program," NASA Contractor Report 187109, Contract NAS3-25556, May 1991. 

2 Goetz, O. K., and Monk, J. C., “Combustion Devices Failures During Space Shuttle Main Engine Development,” Fifth 
International Symposium on Liquid Space Propulsion [CD-ROM], Chattanooga, TN, 2003. 

3 Hulka, J., R., “Scaling of Performance in Liquid Propellant Rocket Engine Combustors,” 20 th Anniversary Meeting of the 
Northern Branch of the Japan Society of Aeronautical and Space Sciences, Sendai, Japan, March 7-9, 2007. 

4 Tucker, K., West, J., Williams, R., Rocker, M., Canabal, F., Robles, B., Garcia, R., and Chenoweth, J., “Using CFD as a Rocket 
Injector Design Tool: Recent Progress at Marshall Space Flight Center,” Fifth International Symposium on Liquid Space Propulsion 
[CD-ROM], Chattanooga, TN, 2003. 

5 Tucker, P. K., Rybak, J. A., Hulka, J. R., Jones, G. W., Nesman, T. and West, J. S., “The NASA Constellation University 
Institutes Project: Thrust Chamber Assembly Virtual Institute,” AIAA Paper No. 2006-4524, July, 2006. 

6 West, J., Westra, D., Lin, J. and Tucker, K., “Accuracy Quantification of the Loci-CHEM Code for Chamber Wall Heat 
Fluxes in a G0 2 /GH 2 Single Element Injector Model Problem,” Third International Workshop on Rocket Combustion Modeling, 
Paris, France, March 13-15, 2006. 

7 Lin, J., West, J. S., Williams, R. W., Tucker, P. K. and Chenoweth, J. D., “CFD Code Validation of Wall Heat Fluxes for a 
G0 2 /GH 2 Single Element Combustor,” AIAA Paper No. 2005-4524, July, 2005. 

8 West, J., Lin, J., Tucker, K. and Chenoweth, J., “Steady State Combustion CFD Analysis of Local Heat Transfer for Liquid 
Oxygen/Gaseous Hydrogen Injectors,” 53 rd JANNAF Propulsion Meeting/2 nd Liquid Propulsion Subcommittee Meeting, 
December 5-8, 2005. 

9 Pal, S., Marshall, W., Woodward, R. and Santoro, R., “Wall Heat Flux Measurements for a Uni-element G0 2 /GH 2 Shear 
Coaxial Injector,” Third International Workshop on Rocket Combustion Modeling, Paris, France, March 13-15, 2006. 

10 Jones, G., Protz, C., Bullard, B., Hulka, J., Santoro, R., Pal, S., Woodward, R. and Marshall, W., “Local Heat Transfer 
Measurements in a L0 2 /GH 2 Single Element Combustor,” 53 rd JANNAF Propulsion Meeting/2 nd Liquid Propulsion 
Subcommittee Meeting, December 5-8, 2005. 

n Oberkampf, W. L., Trucano, T. G., “Verification and Validation in Computational Fluid Dynamics” Progress in Aerospace 
Sciences , Vol. 38, No. 3, pp. 209-272. 

12 Oefelein, J. C., “Large eddy simulation of turbulent combustion processes in propulsion and power systems,” Progress in 
Aerospace Sciences , Vol. 42, 2006, pp. 2-37. 

I3 Oefelein, J. C., “Mixing and combustion of cryogenic oxygen-hydrogen shear-coaxial jet flames at supercritical pressure,” 
Combustion Science and Technology , Vol. 178, No. 1-3, 2006, pp. 229-252. 

14 Oefelein, J. C., “Thermophysical characteristics of LOX-H 2 flames at supercritical pressure,” Proceedings of the 
Combustion Institute , Vol. 30, 2005, pp. 2929-2937. 

15 Pantano, C. and Sarkar, S., “A subgrid model for nonlinear functions of a scalar,” Physics of Fluids, Vol. 13, No. 12, 2001 
pp. 3803-3819. 

16 Mellado, J. P., Sarkar, S., and Pantano, C., “Reconstruction subgrid models for nonpremixed combustion,” Physics of 
Fluids , Vol. 15, No. 1 1, 2003, pp. 3280-3307. 


12 

American Institute of Aeronautics and Astronautics 



17 Conaire, 6. M., Curran, H. J., Simmie, J. M., Pitz, W. J. and Westbrook, C. K., “A comprehensive modeling study of 
hydrogen oxidation,” International Journal of Chemical Kinetics , Vol. 36, 2004, pp. 603-622. 

18 Menon, S. and Patel, N., “Subgrid Combustion Modeling for LES of Spray Combustion in Large-Scale Combustors,” AIAA 
Journal , Vol. 46, No. 4, 2006, pp. 709-723. 

19 E1-Asrag, H. and Menon, S. “Large Eddy Simulation of Bluff-Body Stabilized Swirling Non-Premixed Flames,” 
Proceedings of the Combustion Institute , Vol. 31, 2007, pp. 1747-1754. 

20 Genin, F., Fryxell, B. and Menon, S., “Hybrid Large-Eddy Simulation of Detonation in Reactive Mixtures,” Proceedings of 
the 20th International Conference on Detonations , Explosions and Shock Waves , Montreal CA, August 1-4, 2005. 

21 Meng, H., and Yang, V., “A Unified Treatment of General Fluid Thermodynamics and Its Application to a Preconditioning 
Scheme ,”/. Comput . Phys., Vol. 189, 2003, pp. 277-304. 

22 Meng, H., Hsiao, G.C., Yang, V., and Shuen, J. S., “Transport and Dynamics of Liquid Oxygen Droplets in Supercritical 
Hydrogen Streams,”/. Fluid Mech ., Vol. 527, 2005, pp. 1 15-139. 

23 Zong, N., Meng, H., Hsieh, S.Y., and Yang, V., “A Numerical Study of Cryogenic Fluid Injection and Mixing under 
Supercritical Conditions,” Phys. Fluids , Vol. 16, 2004, pp. 4248-4261. 

24 Yang V., “Modeling of Supercritical Vaporization, Mixing and Combustion Processes in Liquid-Fueled Propulsion 
Systems. Proc. Combust. Inst ., Vol. 28, 2000, pp. 925-942. 

25 Ribert, G., Zong, N., and Yang, V., “Counterflow Diffusion Flames of General Fluids: Oxygen/Hydrogen Mixtures,” AIAA 
Paper No. 2007-1427, 2007. 

26 Zong, N., Ribert, G., and Yang, V., “Supercritical Combustion of Liquid Oxygen (LOX) and Methane Stabilized by a 
Splitter Plate” AIAA Paper No. 2007-0575, 2007. 

27 Hsieh, S.Y., and Yang, V., “A Preconditioned Flux-Differencing Scheme for Chemically Reacting Flows at All Mach 
Numbers,” Int. J. Comput. Fluid Dyn. , Vol. 8, 1997, pp. 31-49. 

28 Zong, N., and Yang, V., “A efficient preconditioning scheme for real fluid mixtures using primitive pressure-temperature 
variables,” submitted to Int. J. Comput. Fluid Dyn., 2007. 

29 Venkateswaran, S., Merkle, C. L., Zeng, X and Li, D., “Influence of Large-Scale Pressure Changes on Pre-Conditioned 
Solutions at Low Speeds,” AIAA Journal , Vol. 42, No. 12, 2004, pp 2490-2498. 

30 Li, D., Xia, G. and Merkle, C. L., “Analysis of Real Fluid Flows in Converging Diverging Nozzles,” AIAA Paper 2003- 
4132, 33rd AIAA Fluid Dynamics Conference and Exhibit, Orlando, Florida, June 23-26, 2003. 

31 Menter, F. R., “Two-Equation Eddy- Viscosity Turbulence Models for Engineering Applications,” AIAA Journal , Vol. 32, 
No. 8, 1994, pp. 1598-1605. 

32 Wilcox, D. C., “Reassessment of the Scale-Determining Equation for Advanced Turbulence Models,” AIAA Journal , Vol. 
26, No. 11, 1988, pp. 1299-1310. 

33 Evans. J. S. and Schexnayder, C. J., “Influence of Chemical Kinetics and Unmixedness on Burning in Supersonic Hydrogen 
Flames,” AIAA Journal , Vol 18, No 2, 1980, pp. 188-193. 

34 Luke, E. A., Tong, X.-L. and Cinnella, P. “Numerical Simulations of Fluids with a General Equation of State”, AIAA 2006- 
1295 1 44th Aerospace Sciences Meeting, January 9-12, 2006, Reno, NV. 

35 Luke, E. A., Tong, X.-L., Wu, J. and Cinella, P., “CHEM 2: A Finite-Rate Viscous Chemistry Solver - The User Guide,” 
MSSU-COE-ERC-04-07, Mississippi State University, September, 2004. 

36 Luke, E. A., “A Rule-Based Specification System for Computational Fluid Dynamics,” Ph.D. Dissertation, Mississippi 
State University, Starkville, MS, 1999. 

37 Gridgen, Software Package, Version 15, Pointwise, Inc., Bedford, Texas, 2005. 


13 

American Institute of Aeronautics and Astronautics 



Table 1, SRL Summary for injector element simulations. 


Element 

Description 

f 

r 


1 

Single Element Concentric Shear Coax G02/GH2 (hot propellants) 

2 

3 

1/1 

2 

Single Element Concentric Shear Coax L02/GH2 

2 

3 

0/3 

3 

Single Element Concentric Swirl Coax L02/GH2 

2 

3 

0/2 
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Table 2: Stream composition, reference properties, and flow characteristics of the oxidizer and fuel streams. 



Composition by Mass 
Composition by Volume 


Reference Properties: 



Pressure, MPa 

5.2 


Temperature, K 

711 

800 

Density, kg/m 3 

26.8 

3.33 

Specific Heat (C p ), J/kg-K 

1110 


Ratio of Specific Heats 

1.34 

1.38 

Dynamic Viscosity, Pa s 

3.62 x 1 0' 5 

1.81 x 10- 5 

Thermal Conductivity, W/m*K 

0.0602 

0.260 

Kinematic Viscosity, m 2 /s 

1.35 x 10' 6 

5.44 x 10' 6 

Thermal Diffusivity, m 2 /s 

2.03 x 10‘ 6 

10.8 x 10' 6 

Sound Speed, m/s 

513 

1470 

Flow Characteristics: 



Bulk Velocity, m/s 

68.0 (inlet) / 154 (exit) 

25.9 (inlet)/ 764 (exit) 

Friction Velocity, m/s 

2.77 (inlet)/ 6.07 (exit) 

1.31 (inlet)/ 35.5 (exit) 

Reynolds Number* 

401,000(inlet) / 604,000 (exit) 

60,000 (inlet) / 169,000 (exit) 

Reynolds Number* 

8170 (inlet)/ 11900 (exit) 

770 (inlet) / 1960 (exit) 


^Based on hydraulic diameter. 
*Based on friction velocity. 
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Figure la. Representative comparison of heat flux predictions from CFD with the corresponding 
experimental data (Element 1: Shear coaxial injector with G0 2 /GH 2 propellants). 



Figure lb. Representative comparison of heat flux predictions from CFD with the corresponding 
experimental data (Element 2: Shear coaxial injector with LQ 2 /GH 2 propellants). 



Figure lc. Representative comparison of heat flux predictions from CFD with the corresponding 
experimental data (Element 3: Swirl coaxial injector with L0 2 /GH 2 propellants). 
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Figure 2. Schematic of the shear coaxial element and combustor test rig. 
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Reference Chamber Pressure: 5.2 MPa 
d = 5.26 mm (Oxidizer Post Inner Diameter) 
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Figure 3. Computational domain and boundary conditions. 
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Figure 2. Normalized instantaneous fuel and oxidizer velocity distributions at 1 dimensionless unit upstream 
of the face plate. Contours represent the axial velocity component, vectors represent the transverse 
components. 
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Figure 5a. Snapshots of near-field temperature distributions. 
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Figure 5b. Snapshots of near-field vorticity magnitude distributions. 
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Figure 5c. Snapshots of near-field hydrogen mass fraction distributions. 
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Figure 5d. Snapshots of near-field oxygen mass fraction distributions. 
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Figure 5e. Snapshots of near-field OH radical mass fraction distributions. 
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Figure 6a. Time-averaged distributions of temperature in the near field. 
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Figure 6b. Time-averaged distributions of vorticity magnitude in the near field. 
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Figure 6c. Time-averaged distributions of turbulent kinetic energy (TKE) in the near field. 
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Figure 6d. Time-averaged distributions of hydrogen mass fraction in the near field. 
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Figure 6e. Time-averaged distributions of oxygen mass fraction in the near field. 
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Figure 6f. Time-averaged distributions of OH radical mass fraction in the near field. 
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Introduction 



♦ NASA plans to modify existing engines to meet Constellation mission 
requirements for Isp 

• Both engines will be modified to provide more thrust 

• A portion of the increased thrust will be derived from higher injector efficiency 
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Introduction 



♦ Injectors on both engines were designed with 0- or 1-D empirical 
tools followed by significant amounts of costly testing 

♦ Injector flows are 

• responsible for TCA performance and thermal environments 

• the result of unsteady, 3-D flows 



MCTA Testing at MSFC 


Injector-generated MCC damage 


Background 



♦ CFD has the potential to replace empiricism in legacy design tools with 
a higher-fidelity, physics-based approach. Accounting for more injector 
physics will lead to: 

• More robust TCA designs via 

Sensitivity analysis for improved conceptual designs 
Unambiguous scaling from subscale test data to full scale design 

• Decreased requirements for large scale testing 

♦ CFD has not been widely used for detailed injector design because of a 
lack of : 

• Model fidelity — geometry & physics 

• Simulation robustness — slow solution process 

• Demonstrated solution accuracy — lack of methodical validation 


Background — Validation Philosophy 
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Background — Previous Validation Efforts 


♦ MSFC has conducted significant single element injector validation work 
over the past 4 years at the single element injector model problem level 

• Coaxial element injectors — shear and swirl 

• 0 2 (GOX & LOX) and GH 2 propellants — at ambient and hot temperatures 

• P r from 300-1200 psia 


Fuel ► 1 set(1 TC & 1 HF) 
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chamber wall temperature profile 



CFD simulated temperature field 



wall heat flux comparison 
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Background — Previous Validation Efforts 



Element 1: Shear coaxial injector 
with G0 2 /GH 2 propellants 


Element 2: Shear coaxial injector 
with L0 2 /GH 2 propellants 


Element 3: Swirl coaxial injector 
with L0 2 /GH 2 propellants 


♦ CFD Model 

• All simulations accomplished with steady, axisymmetric RANS models 

• Comparison to wall heat flux data is fair for Element 1 — under-prediction in the 
near-injector region, over-prediction down stream 

• Model to data comparison degrades significantly with increasing physical and 
geometrical complexity — very significant under-prediction in the near-injector region 

♦ Experiment 

• Heat flux measurement is being used as an indirect indication of mixing 

• There is no guidance for model development in terms of flow field description 
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Background — Current & Future Validation Efforts 


1. Investigate Unit Physics Problems 



Shear 



V//////A 


k k k k 





Heat Transfer 
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Shear, Heat Transfer, 
and Mass Transfer 


2. Continue systematic investigation/validation with Injector Model Problems using 
higher fidelity methods 

• What geometry & physics fidelity are required to match the data? 

• What computing resources are required? 

• What can be achieved in the design/production environment? 
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Scope of the Overall Effort 



♦ Methodologies 

• Sandia (Oefelein) 

• PSU (Yang) 

• Georgia Tech (Menon) 

• Purdue (Merkle) decreased TKE — approaches HRLES 

• MSFC (Tucker) 

♦ Wall Heat Flux Simulations at Injector Model Problem Validation Level 

1. Element 1: Shear coaxial injector with G02/GH 2 propellants 

2. Element 2: Shear coaxial injector with L0 2 /GH 2 propellants 

3. Element 3: Swirl coaxial injector with L0 2 /GH 2 propellants 


10 


Objectives of the Current Paper 


Status initial results for the Element 1 by 
comparing/contrasting simulation results in region just 
downstream of the injector post tip-subset of overall 
calculation 

1 . Extract initial flow field observations 

2. Establish consistencies/inconsistencies in the calculations in terms of: 

a. Model Assembly 

b. Model Implementation 

c. Model Execution 

3. Comment on impact on results of the consistencies/inconsistencies 

4. Provide “best practices” framework to: 

a. Guide future calculations 

b. Facilitate quantitative comparisons that yield clear conclusions to advance the 
validation process and the ability to support injector design issues 


Results— Initial Flow Field Observations 



♦ General flame characteristics 

• Flame holding accomplished by reverse flow in post tip wake 

• Reactions occur along shear layer until 0 2 is consumed 

♦ Unsteady solutions dominated by large vortical motions originating in the post tip 
region and amplified downstream causing significant differences between the time 
averaged (TA) instantaneous solutions and the steady state (SS) solutions: 

• Higher vorticity levels in the TA results indicate increased convective mixing 
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• Propellants in TA results are consumed more quickly than SS results 
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Results — Initial Flow Field Observations (Con’t) 


• Larger region of higher TKE levels in the TA results versus the SS results 
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• Flame 

- Shorter flame length for TA versus SS 

- TA flame is discontinuous versus continuous for SS 

- TA flame is more diffuse and has a lower peak temperature than SS 

- Much hotter gases radially outboard of the injector element for TA versus SS 
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Results-Simulation Consistencies/Inconsistencies 



♦ Model Assembly 

• Boundary condition implementation consistent across different 
models/methodologies 
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• Chemistry model use varied — note OH contours 
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• Turbulence models (for RANS vs URANS) 

• Laminar properties 
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Results-Simulation Consistencies/Inconsistencies 


♦ Model Implementation 

• 2D vs 3D flow field structures are different 

Sandia/3D PSU/2D Purdue/2D 



x (mm) 


• Grids — each technique has its own grid resolution requirements 
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Results-Simulation Consistencies/Inconsistencies 


♦ Model Implementation & Execution 

• Flow field initialization must be consistent to facilitate meaningful comparison of 
instantaneous results 

• Sufficient flow thru times must be achieved to wash out initial condition so that time- 
averaged solutions can be compared to heat flux data 

• PSU Initialization — steady RANS solution interpolated onto LES grid, then run for 0.5 flow 
thru times 

• Purdue Initialization — chamber initialized with quiescient 1500K N 2 , then run for 2.5 flow 
thru times 

PSU instantaneous Purdue instantaneous 


x (mm) 

Purdue time-averaged 


MSFC Steady State 
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Best Practices & Future Work 


♦ Best Practices Identified 

• Model Assembly 

Boundary conditions 
Chemistry mechanisms 
Laminar properties 

• Model Implementation 

Domain — 2D/3D 

Methodology-appropriate grid resolution 
Flow field initialization 

• Model Execution 

Sufficient run time 


♦ Future Work 

• Implement and exercise best practices 

• Repeat initial calculation for Element 1 (Shear coaxial injector with G0 2 /GH 2 propellants) 

• Compare instantaneous results at consistent times in post tip region 

• Compare time-averaged results to experimental heat flux 

• Element 2 (Shear coaxial injector with L0 2 /GH 2 propellants) 

• Element 3 (Swirl coaxial injector with L0 2 /GH 2 propellants) 
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Summary 


♦ Models of Element 1 (G0 2 /GH 2 shear coaxial injector) have been developed using 5 
methodologies ranging from steady RANS to LES 

♦ First solutions have been produced and results compared qualitatively in the critical 
flow field region just downstream of the injector post tip 

♦ In general, the time-averaged results of the unsteady solutions, as compared to the 
steady RANS solutions, exhibit 

• Higher levels of vorticity in the near injector region 

• Faster propellant consumption 

• Flame discontinuity just downstream of the post tip 

• More diffuse flames with lower maximum temperatures 

• Higher temperatures radially outboard of the injector — could be a key in correcting the 
steady RANS under-prediction of the near-injector wall heat flux 

♦ Inconsistencies in model assembly, implementation and execution preclude 
quantitative comparisons at this point 

♦ Best practices in these areas have been identified 

♦ Future work 

• Near term-implementation of best practices 

to continue Element 1 calculations 

to facilitate quantitative comparisons and generate actionable conclusions 

• Longer term — models of Elements 2 and 3 
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